function [B_out] = gen_B(eta, QQ, data,dates0, Ebar,indUnc)

% This is CI Varing B
Vs   = data(:,1);
Gold = data(:,2);
%T = size(Vs,1);

% Parameter
NQ = size(QQ,3);

OMEGA= cov(eta);
BB0=chol(OMEGA,'lower'); %Generate lower cholesy factorization of covariance matrix eta

% dates1      = 200712; %2007:12
% dates2      = 200906; %2009:06
  
  dates3      = 201107; %2011:07
  dates4      = 201108; %2011:08
  
  dates5      = 198710; %1987/10
  dates20     = 197910; %1979:10
% dates7012   = 197012;
  dateslehman = 200809; %2008/09;

% %GFC
% index_07 = find(abs(dates0 - dates1) < 10^(-5)):find(abs(dates0 - dates2) < 10^(-5));

%87 crash
  index_87 =   [find(abs(dates0 - dates5) < 10^(-5))];

% % 1970
% index_70 = [find(abs(dates0 -dates7012)<10^(-6))];

% Lehman Crash
index_lehman = [find(abs(dates0 -dateslehman)<10^(-6)):find(abs(dates0 -dateslehman)<10^(-6))];


%ceiling crisis
index_crsis = [find(dates0 == dates3),find(abs(dates0 -dates4)<10^(-6))];
% 1979
index_79= [find(abs(dates0 -dates20)<10^(-6)):find(abs(dates0 -dates20)<10^(-6))];


% Solve B
parfor iii = 1:NQ
    
    % rotate and normalize B
    B_initial=BB0*QQ(:,:,iii);
    B_est=B_initial*diag(sign(diag(B_initial))); %normalize so that self-irfs are positive
    ehat = eta(:,:)*inv(B_est)';
    
    
    % Compute Correlations
   %phi1      = corr(Vs, ehat(:,1));
    phi2      = corr(Vs, ehat(:,indUnc));
   %phi1_gold = corr(Gold, ehat(:,1));
    phi2_gold = corr(Gold, ehat(:,indUnc));
    
    
    %============ Event constraint ==============
    % gbar1, gbar2, gbar3
    %ind_ef_econ(iii) =   sum(ehat(index_07,2)) <= 0 &(ehat(index_87,indUnc) >= Ebar(1))  &  (ehat(index_lehman,3) >= Ebar(2));
     ind_ef_econ(iii) =                               (ehat(index_87,indUnc) >= Ebar(1))  &  (ehat(index_lehman,indUnc) >= Ebar(2));
    
    % gbar4, gbar5 and gbar 6
    %ind_ef_event_sign(iii) =  (ehat(index_79,1) >= 0 & ehat(index_79,indUnc) >= 0)& (min(ehat(index_crsis,1)) >= 0&min(ehat(index_crsis,indUnc)) >= 0);
     ind_ef_event_sign(iii) =  (                        ehat(index_79,indUnc) >= 0)& (                              min(ehat(index_crsis,indUnc)) >= 0);
    
    % =========== External Variable Constraint ===========
   %ind_ev(iii) =  phi1 <= 0 & phi2 <= 0 & phi1_gold >= 0 & phi2_gold >= 0;
    ind_ev(iii) =              phi2 <= 0                  & phi2_gold >= 0;
    

    ind_ef(iii) = ind_ev(iii).* ind_ef_event_sign(iii).*ind_ef_econ(iii);
    B_store(:,:,iii) = B_est;
    
end

disp('Number of solutions')
disp(sum(ind_ef))

B_out = B_store(:,:,logical(ind_ef));

end